import numpy as np
import cvxpy as cvx
import matplotlib.pyplot as plt
import cupy as cp


class Pattern:
    # 通用参数
    Re = 6371e3  # 地球半径
    c = 3e8  # 光速
    G = 6.67e-11  # 引力常量
    M = 5.965e24  # 地球质量

    def __init__(self, antenna_parameter: dict, radar_parameter: dict, mission_parameter: dict, sample_parameter: dict):
        # 天线参数读入
        self.Nx = antenna_parameter['Nx']
        self.Ny = antenna_parameter['Ny']
        self.Wxt = antenna_parameter['Wxt']
        self.Wyt = antenna_parameter['Wyt']
        self.Wxr = antenna_parameter['Wxr']
        self.Wyr = antenna_parameter['Wyr']
        self.dxRatio = antenna_parameter['dxRatio']
        self.dyRatio = antenna_parameter['dyRatio']
        self.A = antenna_parameter['A']
        self.A_r = np.deg2rad(self.A)

        # 雷达参数读入
        self.H = radar_parameter['H']
        self.f0 = radar_parameter['f0']
        self.PRF = radar_parameter['PRF']
        self.Br = radar_parameter['Br']
        self.N = radar_parameter['N']

        # 任务参数读入
        self.alpha_c = mission_parameter['alphaC']
        self.theta_sq_c = mission_parameter['thetaSqC']
        self.alpha_scan_min = mission_parameter['alphaScanMin']
        self.alpha_scan_max = mission_parameter['alphaScanMax']
        self.alpha_c_r = np.deg2rad(self.alpha_c)
        self.theta_sq_c_r = np.deg2rad(self.theta_sq_c)
        self.alpha_scan_max_r = np.deg2rad(self.alpha_scan_max)
        self.alpha_scan_min_r = np.deg2rad(self.alpha_scan_min)

        # 采样手动设置
        self.R_sample_rate = sample_parameter['RSampleRate']
        self.f_sample_rate = sample_parameter['fSampleRate']

        # 雷达参数计算
        self.PRF_r = self.PRF * self.N
        self.pr = self.c / 2 / self.Br
        self.lamda = self.c / self.f0  # 波长
        self.v = (self.G * self.M / (self.H + self.Re)) ** 0.5  # 雷达平台速度
        self.f_eta_c = 2 * self.v * np.sin(self.theta_sq_c_r) / self.lamda  # 多普勒中心频率

        # 天线参数
        self.dx = self.dxRatio * self.lamda  # 天线x单元间隔
        self.dy = self.dyRatio * self.lamda  # 天线y单元间隔

        # 最大扫描角度参数计算
        self.alpha_ant_min = self.A - 90  # 雷达最小可视角度
        self.alpha_ant_max = self.A + 90  # 雷达最大可视角度
        self.alpha_ant_min_r = np.deg2rad(self.alpha_ant_min)
        self.alpha_ant_max_r = np.deg2rad(self.alpha_ant_max)
        self.theta_earth_r = np.arcsin(self.Re / (self.Re + self.H))  # 地球最大扫描角度
        self.theta_earth = np.rad2deg(self.theta_earth_r)

        # 主信号区
        self.theta_min_r = self.cal_theta(self.alpha_scan_min_r, self.theta_sq_c_r)
        self.theta_max_r = self.cal_theta(self.alpha_scan_max_r, self.theta_sq_c_r)
        self.theta_min = np.rad2deg(self.theta_min_r)
        self.theta_max = np.rad2deg(self.theta_max_r)

        # 雷达姿态
        self.A = np.arctan(np.tan(self.alpha_c_r) / np.cos(self.theta_sq_c_r))
        self.theta_yaw = np.arctan(np.sin(self.alpha_c_r) * np.tan(self.theta_sq_c_r))
        self.theta_pit = np.arcsin(np.cos(self.alpha_c_r) * np.sin(self.theta_sq_c_r))

        # 多普勒
        self.theta_bw_r = 0.886 / self.dyRatio / self.Ny  # 波束宽度
        self.theta_bw = np.rad2deg(self.theta_bw_r)
        self.Bd = 2 * self.v * np.cos(self.theta_sq_c_r) * self.theta_bw_r / self.lamda  # 多普勒带宽

        self.R_min = self.theta2r_eta_cpu(self.theta_min_r)[0]  # 最近斜距
        self.R_max = self.theta2r_eta_cpu(self.theta_max_r)[0]  # 最远斜距

        # 计算模糊距离
        self.R0 = np.linspace(self.R_min, self.R_max, self.R_sample_rate)

    def theta2r_eta_cpu(self, theta):  # 由星下点角计算斜距、入射角
        eta = np.arcsin(np.round(np.sin(theta) * (self.Re + self.H) / self.Re, 12))
        beta = eta - theta
        r = np.sin(beta) / np.sin(theta) * self.Re
        if type(r) == np.float64 and r == 0:
            r = self.H
        elif type(r) == np.array:
            r[np.where(beta == 0)] = self.H
        return r, eta

    def steer(self, alpha, theta_sq):
        shape1 = list(alpha.shape)
        shape1.extend([self.Nx, 1])
        shape2 = list(theta_sq.shape)
        shape2.extend([self.Ny, 1])

        ex = cp.zeros(tuple(shape1), dtype='complex')
        ey = cp.zeros(tuple(shape2), dtype='complex')
        for i in range(self.Nx):
            ex[..., i, 0] = cp.exp(2j * cp.pi / self.lamda * i * self.dx *
                                   (cp.cos(theta_sq) * cp.sin(alpha - self.A_r) -
                                    cp.cos(self.theta_sq_c_r) * cp.sin(self.alpha_c_r - self.A_r)))
        for i in range(self.Ny):
            ey[..., i, 0] = cp.exp(2j * cp.pi / self.lamda * i * self.dy *
                                   (cp.sin(theta_sq) - cp.sin(self.theta_sq_c_r)))

        trans = [i for i in range(len(ex.shape))]
        trans.insert(-1, trans.pop())

        e = ex @ cp.transpose(ey, tuple(trans))
        shape = list(alpha.shape)
        shape.extend([self.Nx * self.Ny, 1])
        return e.reshape(tuple(shape))

    def steer_cpu(self, alpha, theta_sq):
        shape1 = list(alpha.shape)
        shape1.extend([self.Nx, 1])
        shape2 = list(theta_sq.shape)
        shape2.extend([self.Ny, 1])

        ex = np.zeros(tuple(shape1), dtype='complex')
        ey = np.zeros(tuple(shape2), dtype='complex')
        for i in range(self.Nx):
            ex[..., i, 0] = np.exp(2j * np.pi / self.lamda * i * self.dx *
                                   (np.cos(theta_sq) * np.sin(alpha - self.A_r) -
                                    np.cos(self.theta_sq_c_r) * np.sin(self.alpha_c_r - self.A_r)))
        for i in range(self.Ny):
            ey[..., i, 0] = np.exp(2j * np.pi / self.lamda * i * self.dy *
                                   (np.sin(theta_sq) - np.sin(self.theta_sq_c_r)))

        trans = [i for i in range(len(ex.shape))]
        trans.insert(-1, trans.pop())

        e = ex @ np.transpose(ey, tuple(trans))
        shape = list(alpha.shape)
        shape.extend([self.Nx * self.Ny, 1])
        return e.reshape(tuple(shape))

    def steer_2d(self, alpha, theta_sq):
        def phase(al, th):
            phx = ((cp.cos(self.A) * cp.cos(self.theta_yaw) +
                    cp.sin(self.A) * cp.sin(self.theta_yaw) * cp.sin(self.theta_pit)) * cp.sin(al) * cp.cos(th) +
                   (cp.cos(self.A) * cp.sin(self.theta_yaw) -
                    cp.sin(self.A) * cp.cos(self.theta_yaw) * cp.sin(self.theta_pit)) * cp.sin(th) -
                   cp.sin(self.A) * cp.cos(self.theta_pit) * cp.cos(al) * cp.cos(th)) * \
                  self.dx * 2j * cp.pi / self.lamda

            phy = (-cp.cos(self.theta_pit) * cp.sin(self.theta_yaw) * cp.sin(al) * cp.cos(th) +
                   cp.cos(self.theta_pit) * cp.cos(self.theta_yaw) * cp.sin(th) -
                   cp.sin(self.theta_pit) * cp.cos(al) * cp.cos(th)) * self.dy * 2j * cp.pi / self.lamda

            phz = ((cp.sin(self.A) * cp.cos(self.theta_yaw) -
                    cp.cos(self.A) * cp.sin(self.theta_yaw) * cp.sin(self.theta_pit)) * cp.sin(al) * cp.cos(th) +
                   (cp.sin(self.A) * cp.sin(self.theta_yaw) +
                    cp.cos(self.A) * cp.cos(self.theta_yaw) * cp.sin(self.theta_pit)) * cp.sin(th) +
                   cp.cos(self.A) * cp.cos(self.theta_pit) * cp.cos(al) * cp.cos(th))
            return phx, phy, phz

        phase_x, phase_y, phase_z = phase(alpha, theta_sq)

        shape_x = list(phase_x.shape)
        shape_y = list(phase_y.shape)
        shape_x.extend([self.Nx, 1])
        shape_y.extend([self.Ny, 1])
        sx = cp.zeros(tuple(shape_x), dtype='complex')
        sy = cp.zeros(tuple(shape_y), dtype='complex')
        p = cp.zeros(phase_x.shape, dtype='complex')
        for i in range(4):
            p += cp.exp(i * phase_y)
        for i in range(self.Nx):
            sx[..., i, 0] = cp.exp(i * phase_x)
        for i in range(self.Ny):
            sy[..., i, 0] = p * cp.exp(4 * i * phase_y)

        # loc = np.where(phase_z < 0)
        # sx[loc[0], loc[1], ...] = cp.nan + cp.nan * 1j
        # sy[loc[0], loc[1], ...] = cp.nan + cp.nan * 1j

        return sx, sy

    def vector_2d(self, alpha, theta_sq):
        sx, sy = self.steer_2d(alpha, theta_sq)
        v = self.vector(sx, sy)
        shape = list(alpha.shape)
        shape.extend([self.Nx * self.Ny, 1])
        v = v.reshape(shape)
        return v

    def e2m(self, s):
        m = np.zeros([s.shape[0], self.Nx * self.Ny, self.Nx * self.Ny], dtype='complex')
        t = s.copy()
        for i in range(self.Nx):
            m[:, self.Ny * i: self.Ny * (i + 1), :] = t
            if i < self.Nx - 1:
                t[:, :, self.Ny:] = t[:, :, 0: (self.Nx - 1) * self.Ny]
                t[:, :, 0: self.Ny] = np.transpose(np.conj(s[:, :, self.Ny * (i + 1): self.Ny * (i + 2)]), (0, 2, 1))
        return m

    def plot_pattern(self, wmn, if3d=False):
        alpha = np.linspace(self.alpha_ant_min, self.alpha_ant_max, 361)
        theta_sq = np.linspace(-90, 90, 361)
        alpha_r = np.deg2rad(alpha)
        theta_sq_r = np.deg2rad(theta_sq)
        alpha, theta_sq = np.meshgrid(alpha, theta_sq)
        alpha_r, theta_sq_r = np.meshgrid(alpha_r, theta_sq_r)
        e = self.steer_cpu(alpha_r, theta_sq_r)
        s = wmn.T @ e
        s = s.reshape(alpha.shape)
        s_db = abs(s)
        s_db = s_db / np.max(s_db)
        s_db = 20 * np.log10(s_db)
        # s_db_raw = np.copy(s_db)
        s_db[np.where(s_db < -50)] = -50

        if if3d:
            ax = plt.axes(projection='3d')
            ax.plot_surface(alpha.get(), theta_sq.get(), s_db.get(), rstride=5, cstride=5)
            ax.set_xlabel(r'$\alpha$')
            ax.set_ylabel(r'$\theta_{sq}$')
            ax.set_zlabel('S\\dB')
            ax.set_zlim([-50, 0])
            plt.figure()

        plt.contourf(alpha, theta_sq, s_db, cmap='jet')
        plt.colorbar()
        plt.xlabel(r'$\alpha$')
        plt.ylabel(r'$\theta_{sq}$')
        plt.figure()
        plt.pcolor(alpha, theta_sq, s_db, cmap='jet', vmax=0, vmin=-50)
        plt.colorbar()
        plt.xlabel(r'$\alpha$')
        plt.ylabel(r'$\theta_{sq}$')

    def plot_pattern_rotate(self, wmn, if3d=False):
        alpha = cp.linspace(self.alpha_ant_min, self.alpha_ant_max, 361)
        theta_sq = cp.linspace(-90, 90, 361)
        alpha_r = cp.deg2rad(alpha)
        theta_sq_r = cp.deg2rad(theta_sq)
        alpha, theta_sq = cp.meshgrid(alpha, theta_sq)
        alpha_r, theta_sq_r = cp.meshgrid(alpha_r, theta_sq_r)
        e = self.vector_2d(alpha_r, theta_sq_r)
        s = wmn.T @ e
        s = s.reshape(alpha.shape)
        s_db = abs(s)
        s_db = cp.asnumpy(s_db)
        s_db = s_db / np.nanmax(s_db)
        s_db = 20 * np.log10(s_db)
        # s_db_raw = np.copy(s_db)
        s_db[np.where(s_db < -50)] = -50

        if if3d:
            ax = plt.axes(projection='3d')
            ax.plot_surface(alpha.get(), theta_sq.get(), s_db, rstride=5, cstride=5)
            ax.set_xlabel(r'$\alpha$')
            ax.set_ylabel(r'$\theta_{sq}$')
            ax.set_zlabel('S\\dB')
            ax.set_zlim([-50, 0])
            plt.figure()

        plt.contourf(alpha.get(), theta_sq.get(), s_db, cmap='jet')
        plt.colorbar()
        plt.xlabel(r'$\alpha$')
        plt.ylabel(r'$\theta_{sq}$')
        plt.figure()
        plt.pcolor(alpha.get(), theta_sq.get(), s_db, cmap='jet', vmax=0, vmin=-50)
        plt.colorbar()
        plt.xlabel(r'$\alpha$')
        plt.ylabel(r'$\theta_{sq}$')

    def plot_rar_r_separated(self, ss, sa, wx, wy, if_plot=False):
        wmn = wx @ wy.T
        wmn = wmn.reshape(self.Nx * self.Ny, 1)
        ss = wmn.T @ ss @ np.conj(wmn)
        sa = wmn.T @ sa @ np.conj(wmn)
        rar = (abs(sa) / abs(ss)).reshape(self.R_sample_rate, )
        rar_db = 10 * np.log10(rar)
        if if_plot:
            plt.plot(self.R0, rar_db)
            plt.xlabel('distance')
            plt.ylabel('RASR/dB')
        return rar_db

    def plot_rar_r(self, ss, sa, wmn, if_plot=False):
        wmn = wmn.reshape(self.Nx * self.Ny, 1)
        ss = wmn.T @ ss @ np.conj(wmn)
        sa = wmn.T @ sa @ np.conj(wmn)
        rar = (abs(sa) / abs(ss)).reshape(self.R_sample_rate, )
        rar_db = 10 * np.log10(rar)
        if if_plot:
            plt.plot(self.R0, rar_db)
            plt.xlabel('distance')
            plt.ylabel('RASR/dB')
        return rar_db

    @staticmethod
    def cal_alpha_all(theta_sq, theta):
        theta_sq, theta = cp.meshgrid(theta_sq, theta)
        return cp.arccos(cp.cos(theta) / cp.cos(theta_sq)), theta_sq, theta

    @staticmethod
    def cal_theta(alpha, theta_sq):  # input:radian output:radian
        return np.arccos(np.cos(alpha) * np.cos(theta_sq))

    @staticmethod
    def vector(sx, sy):
        sy = cp.transpose(sy, (0, 1, 3, 2))
        v = sx @ sy
        return v

    def doppler2theta_sq(self, fd):
        return np.arcsin(self.lamda * fd / 2 / self.v)


class ConvexProblem:
    def __init__(self, antenna_parameter, pattern_info, va_info):
        # 模糊矩阵
        self.S0 = va_info['S0']
        self.Sa = va_info['Sa']

        self.Sa = self.c2r(self.Sa)
        self.S0 = self.c2r(self.S0)

        # 导向矢量
        # 主瓣
        self.e_mu = self.e2e(pattern_info['Mu'])
        self.e_md = self.e2e(pattern_info['Md'])
        self.e_c = self.e2e(pattern_info['C'])

        # 副瓣
        self.e_s = self.e2e(pattern_info['S'])

        # 天线参数
        self.Nx = antenna_parameter['Nx']
        self.Ny = antenna_parameter['Ny']

    @staticmethod
    def c2r(cm):
        return np.block([[cm.real, -cm.imag], [cm.imag, cm.real]])

    @staticmethod
    def e2e(e):
        return np.concatenate((e, 1j * e), axis=1)

    def problem_solve(self):
        w = cvx.Variable(2 * self.Nx * self.Ny)
        cvx.settings.EIGVAL_TOL = 1e-4
        obj = cvx.Minimize((1/2) * cvx.quad_form(w, self.Sa))
        constrains = []
        constrains += [cvx.abs(self.e_s[i, ...].T @ w) <= 10 ** (20 / 20) for i in range(self.e_s.shape[0])]
        # constrains += [cvx.abs(self.e_md[i, ...].T @ w) >= 10 ** (-3 / 20) for i in range(self.e_md.shape[0])]
        # constrains += [cvx.abs(self.e_mu[i, ...].T @ w) <= 1 for i in range(self.e_mu.shape[0])]
        constrains += [self.e_c[0, ...].T @ w == 1]
        prob = cvx.Problem(obj, constrains)
        prob.solve(solver=cvx.SCS)
        print(prob.status)
        print(w.value)
        return w.value

    def run(self):
        return self.problem_solve()


if __name__ == '__main__':
    AntennaParameter = {
        'Nx': 22,
        'Ny': 16,
        'Wxt': cp.ones((22, 1), dtype='complex'),
        'Wyt': cp.ones((16, 1), dtype='complex'),
        'Wxr': cp.ones((22, 1), dtype='complex'),
        'Wyr': cp.ones((16, 1), dtype='complex'),
        'dxRatio': 0.65,
        'dyRatio': 0.64,
        'A': 30.5
    }

    RadarParameter = {
        'H': 607e3,
        'f0': 1.26e9,
        'PRF': 3250,
        'Br': 80e6,
        'N': 1  # 通道数
    }

    MissionParameter = {
        'alphaC': 26,
        'thetaSqC': 15,
        'alphaScanMin': 25.06,
        'alphaScanMax': 27.33
    }

    SampleParameter = {
        'RSampleRate': 50,
        'fSampleRate': 5
    }

    pattern = Pattern(AntennaParameter, RadarParameter, MissionParameter, SampleParameter)

    interval1 = np.pi / 36
    interval2 = np.pi / 50
    # interval2 = 0
    theta_sq_max = pattern.theta_earth
    theta_sq_min = -pattern.theta_earth
    theta_scan_min = pattern.theta_min
    theta_scan_max = pattern.theta_max
    theta_sq_min_r = np.deg2rad(theta_sq_min)
    theta_sq_max_r = np.deg2rad(theta_sq_max)
    theta_scan_min_r = pattern.theta_min_r
    theta_scan_max_r = pattern.theta_max_r

    f_eta_c = pattern.f_eta_c
    Bd = pattern.Bd
    f_eta_max = f_eta_c + Bd / 2
    f_eta_min = f_eta_c - Bd / 2

    theta_sq_scan_max_r = pattern.doppler2theta_sq(f_eta_max)
    theta_sq_scan_min_r = pattern.doppler2theta_sq(f_eta_min)
    theta_sq_scan_min = np.rad2deg(theta_sq_scan_min_r)
    theta_sq_scan_max = np.rad2deg(theta_sq_scan_max_r)

    alpha_ant_min = pattern.alpha_ant_min
    alpha_ant_max = pattern.alpha_ant_max
    alpha_ant_min = -pattern.theta_earth
    alpha_ant_max = pattern.theta_earth
    alpha_ant_min_r = np.deg2rad(alpha_ant_min)
    alpha_ant_max_r = np.deg2rad(alpha_ant_max)

    # 主信号区域
    theta_sq_main_down = cp.linspace(theta_sq_scan_min_r, theta_sq_scan_max_r, 4)
    theta_main_down = cp.linspace(theta_scan_min_r, theta_scan_max_r, 4)
    Alpha_md, Theta_sq_md, _ = pattern.cal_alpha_all(theta_sq_main_down, theta_main_down)
    steer_main_down = pattern.vector_2d(Alpha_md, Theta_sq_md).reshape(-1, pattern.Nx * pattern.Ny, 1)

    theta_sq_main_up = cp.linspace(theta_sq_scan_min_r - interval2, theta_sq_scan_max_r + interval2, 5)
    theta_main_up = cp.linspace(theta_scan_min_r - interval1, theta_scan_max_r + interval1, 6)
    Alpha_mu, Theta_sq_mu, _ = pattern.cal_alpha_all(theta_sq_main_up, theta_main_up)
    steer_main_up = pattern.vector_2d(Alpha_mu, Theta_sq_mu).reshape(-1, pattern.Nx * pattern.Ny, 1)

    # 旁瓣区域
    # 上
    theta_sq_up = cp.linspace(theta_sq_scan_max_r + interval2, theta_sq_max_r,
                              int((theta_sq_max_r - theta_sq_scan_max_r - interval2) / (np.pi / 60)))
    alpha_up = cp.linspace(alpha_ant_min_r, alpha_ant_max_r,
                           int((pattern.alpha_ant_max_r - pattern.alpha_ant_min_r) / (np.pi / 90)))
    Alpha_up, Theta_sq_up = cp.meshgrid(alpha_up, theta_sq_up)
    steer_up = pattern.vector_2d(Alpha_up, Theta_sq_up).reshape(-1, pattern.Nx * pattern.Ny, 1)

    # 下
    theta_sq_down = cp.linspace(theta_sq_min_r, theta_sq_scan_min_r - interval2,
                                int((theta_sq_scan_min_r - interval2 - theta_sq_min_r) / (np.pi / 60)))
    alpha_down = cp.linspace(alpha_ant_min_r, alpha_ant_max_r,
                             int((pattern.alpha_ant_max_r - pattern.alpha_ant_min_r) / (np.pi / 90)))
    Alpha_down, Theta_sq_down = cp.meshgrid(alpha_down, theta_sq_down)
    steer_down = pattern.vector_2d(Alpha_down, Theta_sq_down).reshape(-1, pattern.Nx * pattern.Ny, 1)

    # 左
    Alpha_left = cp.array([cp.linspace(alpha_ant_min_r, Alpha_mu[0, i],
                                       int((Alpha_mu[0, 0] - alpha_ant_min_r) / (cp.pi / 90)))
                           for i in range(Alpha_mu.shape[1])])
    _, Theta_sq_left = cp.meshgrid(cp.ones(Alpha_left.shape[1]), theta_sq_main_up)
    steer_left = pattern.vector_2d(Alpha_left, Theta_sq_left).reshape(-1, pattern.Nx * pattern.Ny, 1)

    # 右
    Alpha_right = cp.array([cp.linspace(Alpha_mu[-1, i], alpha_ant_max_r,
                                        int((alpha_ant_max_r - Alpha_mu[-1, -1]) / (cp.pi / 90)))
                            for i in range(Alpha_mu.shape[1])])
    _, Theta_sq_right = cp.meshgrid(cp.ones(Alpha_right.shape[1]), theta_sq_main_up)
    steer_right = pattern.vector_2d(Alpha_right, Theta_sq_right).reshape(-1, pattern.Nx * pattern.Ny, 1)

    steer_side = cp.concatenate((steer_up, steer_down, steer_left, steer_right))
    del steer_up
    del steer_down
    del steer_left
    del steer_right

    # 中
    alpha_c_r = cp.array(pattern.alpha_c_r).reshape(1, )
    theta_sq_c_r = cp.array(pattern.theta_sq_c_r).reshape(1, )
    alpha_c_r, theta_sq_c_r = cp.meshgrid(alpha_c_r, theta_sq_c_r)
    steer_c = pattern.vector_2d(alpha_c_r, theta_sq_c_r)
    steer_c = steer_c.reshape(-1, pattern.Nx * pattern.Ny, 1)

    PatternInfo = {
        'Mu': cp.asnumpy(steer_main_up),
        'Md': cp.asnumpy(steer_main_down),
        'S': cp.asnumpy(steer_side),
        'C': cp.asnumpy(steer_c)
    }

    S0, Sa = np.load('S026.npy'), np.load('Sa26.npy')
    S0, Sa = pattern.e2m(S0), pattern.e2m(Sa)
    S0s, Sas = np.sum(S0, axis=0), np.sum(Sa, axis=0)
    S0m = np.max(S0s)
    Sam = np.max(Sas)

    # Sas += 0.2 * Sam * np.eye(pattern.Nx * pattern.Ny, dtype='complex')

    RarInfo = {
        'S0': S0s,
        'Sa': Sas
    }

    Spectral = ConvexProblem(AntennaParameter, PatternInfo, RarInfo)
    w_opt = Spectral.run()

    wmn_opt = w_opt[0:pattern.Nx*pattern.Ny] + 1j * w_opt[pattern.Nx*pattern.Ny:]
    wmn_opt = wmn_opt.reshape(-1, 1)
    plt.figure()

    rar_opt = pattern.plot_rar_r(S0, Sa, wmn_opt, if_plot=True)
    # plt.figure()
    rar_ini = pattern.plot_rar_r_separated(S0, Sa, pattern.Wxr.get(), pattern.Wyr.get(), if_plot=True)
    plt.figure()
    plt.plot(pattern.R0, rar_opt, label='opt')
    plt.plot(pattern.R0, rar_ini, label='ini')
    plt.legend()
    # plt.figure()
    # pattern.plot_rar_theta(S0, Sa, wx_opt, wy_opt, if_plot=True)
    plt.figure()
    pattern.plot_pattern_rotate(cp.asarray(wmn_opt))
    plt.figure()
    pattern.plot_pattern_rotate(cp.ones((pattern.Nx*pattern.Ny, 1)))
    plt.show(block=True)

    print('Program End')